home *** CD-ROM | disk | FTP | other *** search
/ Magnum One / Magnum One (Mid-American Digital) (Disc Manufacturing).iso / d18 / nrpas13.arc / POWELL.PAS < prev    next >
Pascal/Delphi Source File  |  1991-05-01  |  1KB  |  56 lines
  1. PROCEDURE powell(VAR p: glnarray; VAR xi: glnpbynp; n,np: integer;
  2.        ftol: real; VAR iter: integer; VAR fret: real);
  3. (* Programs using routine POWELL must define the types
  4. TYPE
  5.    glnarray = ARRAY [1..n] OF real;
  6.    glnpbynp = ARRAY [1..np,1..np] OF real;
  7. in the main routine. *)
  8. LABEL 1,99;
  9. CONST
  10.    itmax=200;
  11. VAR
  12.    j,ibig,i: integer;
  13.    t,fptt,fp,del: real;
  14.    pt,ptt,xit: glnarray;
  15. BEGIN
  16.    fret := fnc(p);
  17.    FOR j := 1 TO n DO BEGIN
  18.       pt[j] := p[j]
  19.    END;
  20.    iter := 0;
  21. 1: iter := iter+1;
  22.    fp := fret;
  23.    ibig := 0;
  24.    del := 0.0;
  25.    FOR i := 1 TO n DO BEGIN
  26.       FOR j := 1 TO n DO BEGIN
  27.          xit[j] := xi[j,i]
  28.       END;
  29.       fptt := fret;
  30.       linmin(p,xit,n,fret);
  31.       IF (abs(fptt-fret) > del) THEN BEGIN
  32.          del := abs(fptt-fret);
  33.          ibig := i
  34.       END
  35.    END;
  36.    IF (2.0*abs(fp-fret) <= ftol*(abs(fp)+abs(fret))) THEN GOTO 99;
  37.    IF (iter = itmax) THEN BEGIN
  38.       writeln('pause in routine POWELL');
  39.       writeln('too many interations'); readln
  40.    END;
  41.    FOR j := 1 TO n DO BEGIN
  42.       ptt[j] := 2.0*p[j]-pt[j];
  43.       xit[j] := p[j]-pt[j];
  44.       pt[j] := p[j]
  45.    END;
  46.    fptt := fnc(ptt);
  47.    IF (fptt >= fp) THEN GOTO 1;
  48.    t := 2.0*(fp-2.0*fret+fptt)*sqr(fp-fret-del)-del*sqr(fp-fptt);
  49.    IF (t >= 0.0) THEN GOTO 1;
  50.    linmin(p,xit,n,fret);
  51.    FOR j := 1 TO n DO BEGIN
  52.       xi[j,ibig] := xit[j]
  53.    END;
  54.    GOTO 1;
  55. 99:   END;
  56.